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Abstract. Interpolation-based techniques have been widely and successfully ap- 
plied in the verification of hardware and software, e.g., in bounded-model check- 
ing, CEGAR, SMT, etc., whose hardest part is how to synthesize interpolants. 
Various work for discovering interpolants for propositional logic, quantifier-free 
fragments of first-order theories and their combinations have been proposed. 
However, little work focuses on discovering polynomial interpolants in the lit- 
erature. In this paper, we provide an approach for constructing non-linear inter- 
polants based on semidefinite programming, and show how to apply such results 
to the verification of programs by examples. 

Keywords: Craig interpolant, Positivstellensatz Theorem, semidefinite program- 
ming, program verification. 



1 Introduction 

It becomes a grand challenge how to guarantee the correctness of software, as our 
modern life more and more depends on computerized systems. There are lots of ver- 
ification techniques based on either model-checking [1 1, or theorem proving M2I3I . or 
abstract interpretation [4|, or their combination, that have been invented for the verifi- 
cation of hardware and software, like bounded model-checking [5 1, CEGAR [6|, sat- 
isfiability modulo theories (SMT) [7|, etc. The bottleneck of these techniques is scal- 
ability, as many of real software are very complex with different features like compli- 
cated data structures, concurrency, distributed, real-time and hybrid, and so on. While 
interpolation-based techniques provide a powerful mechanism for local and modular 
reasoning, which indeed improves the scalability of these techniques, in which the no- 
tion of Craig interpolants plays a key role. 

Interpolation-based local and modular reasoning was first applied in theorem prov- 
ing due to Nelson and Oppen 0, called Nelson-Oppen method. The basic idea of 
Nelson-Oppen method is to reduce the satisfiability (validity) of a composite theory 
into the ones of its component theories whose satisfiability (validity) have been ob- 
tained. The hardest part of the method, which also determines the efficiency of the 
method, is to construct a formula using the common part of the component theories 
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for a given formula of the composite theory with Craig Interpolation Theorem. In the 
past decade, the Nelson-Oppen method was further extended to SMT which is based 
on DPLL |9| and Craig Interpolation Theorem [ 1 1 for combining different decision 
procedures in order to verify a property of programs with complicated data structures. 
For instance, Z3 ifTTI integrates more than 10 different decision procedures up to now, 
including propositional logic, equality logic with uninterpreted functions, Presburger 
arithmetic, array logic, difference arithmetic, bit vector logic, and so on. 

In recent years, it is noted that interpolation based local and modular reasoning is 
quite helpful for improving the scalability of model-checking, in particular for bounded 
model-checking of systems with finite or infinite states 051121131 . CEGAR 11141151 . etc. 
McMillian first considered how to combine Craig interpolants with bounded model- 
checking to verify infinite state systems. The basic idea of his approach is to generate 
invariant using Craig interpolants so that it can be claimed that an infinite state system 
satisfies a property after k steps in model-checking whenever an invariant is obtained 
which is strong enough to guarantee the property. While in M14I15I16I . how to apply 
the local property of Craig interpolants generated from a counter-example to refine the 
abstract model in order to exclude the spurious counter-example in CEGAR was inves- 
tigated. Meanwhile, in |17|, using interpolation technique to generate a set of atomic 
predicates as the base of machine-learning based verification technique was investigated 
by Wang et al. 

Obviously, synthesizing Craig interpolants is the cornerstone of interpolation based 
techniques. In fact, many approaches have been proposed in the literature. In ifOI . 
McMillian presented a method for deriving Craig interpolants from proofs in the quantifier- 
free theory of linear inequality and uninterpreted function symbols, and based on which 
an interpolating theorem prover was provided. For improving the efficiency of con- 
structing interpolant, McMillian further proposed a method based on lazy abstraction 
for generating interpolants. While, in fl31 . Henzinger et al. proposed a method to syn- 
thesizing Craig interpolants for a theory with arithmetic and pointer expressions, as 
well as call-by-value functions. In lfl8l . Yorsh and Musuvathi presented a combina- 
tion method for generating Craig interpolants for a class of first-order theories. While 
Rybalchenko and Sofronie-Stokkermans |fl9l proposed an approach by reducing the 
synthesis of Craig interpolants of the combined theory of linear arithmetic and uninter- 
preted function symbols to constraint solving. 

However, in the literature, there is little work on how to synthesize non-linear inter- 
polants, except that in [20 1 Kupferschmid and Becker provided a method to construct 
non-linear Craig Interpolant using iSAT, which is a variant of SMT solver based on 
interval arithmetic. 

In this paper we investigate how to construct non-linear interpolants. The idea of 
our approach is as follows: Firstly, we reduce the problem of generating interpolants 
for arbitrary two polynomial formulas to that of generating interpolants for two semi- 
algebraic systems (SASs), which is a conjunction of a set of polynomial equations, 
inequations and inequalities (see the definition later). Then, according to Positivstel- 
lensatz Theorem of real algebraic geometry [21 1, there exists a witness to indicate the 
considered two SASs do not have common real solutions if their conjunction is unsat- 
isfiable. Parrilo in [22 23 1 gave an approach for constructing the witness by applying 



semidefinite programming ll24l . Our algorithm invokes Parrilo's method as a subrou- 
tine. Our purpose is to construct Craig interpolants, so we need to obtain a special wit- 
ness. In general case, we cannot guarantee the existence of the special witness, which 
means that our approach is only sound, but not complete. However, we discuss that if 
the considered two SASs meet Archimedean condition, (e.g., each variable occurring 
in the SASs is bounded, which is a reasonable assumption in practice), our approach is 
not only sound, but also complete. We demonstrate our approach by some examples, in 
particular, we show how to apply the results to program verification by examples. 

The complexity of our approach is polynomial in ub( n+b ^ 2 ) ( n ^ b ), where u is the 
number of polynomial constraints in the considered problem, n is the number of vari- 
ables, and b is the highest degree of polynomials and interpolants. So, the complexity 
of our approach is polynomial in b for a given problem as in which n and u are fixed. 

Structure of the paper: The rest of the paper is organized as follows. By a running 
example, we sketch our approach and show how to apply it to program verification 
in Section [2] Some necessary preliminaries are introduced in Section [3] A sound but 
incomplete algorithm for synthesizing non-linear interpolants in general case is de- 
scribed in Section [4] Section [5] provides a practical algorithm for systems only con- 
taining non-strict inequalities and satisfying Archimedean condition. Section|6]focuses 
on the correctness and complexity analysis of our approach. Our implementation and 
experimental results are briefly reported in Section|7] Section|8]describes more related 
work related to interpolant generation and its application to program verification. Our 
summarizes the paper and discusses future work in Section|9]. 

2 An Overview of Our Approach 

In this section, we sketch our approach and show how to apply our results to pro- 
gram verification by an example. 

1 if (x*x+y*y<l) 

2 / * initial values 

3 while ( x*x+y*y <3){ 

4 x: = x*x+y — 1; 

5 y:=y+x*y +1; 

6 if ( x*x—2*y*y—4>0) 

7 /* unsafe area 

8 error ( ) ; 

9 } 

Code 1.1: example 



Consider the program in Code |l.l| (left). This program tests the initial value of x and y 
at line 1, afterwards executes the while loop with x 2 + y 2 < 3 as the loop condition. The 
body of the while loop contains two assignments and an if statement in sequence. The 



gi = 1 - x 2 — y 2 > 

, 92 = 3 - x 2 - y 2 > 
fx = x 2 + y - 1 - x 1 = 
h = v + x'y + l - y' - 
.93 = x' 2 - 2y' 2 - 4 > 



property we wish to check is that error() procedure will never be executed. Suppose 
there is an execution 1— >3— >4— !>5— !>6— > 8. We can encode such an execution 
by the formulas as in Code (right). Note that in these formulas we use unprimed 
and primed versions of each variable to represent the values of the variable before and 
after updating respectively. Obviously, the execution is infeasible iff the conjunction of 
these formulas is unsatisfiable. Let ^>=<7i>0A/i=0A/2 = (Qand if, = g 3 > 0. 
To show <j> A tp is unsatisfiable, we need to construct an interpolant 9 for tfi and ift, i.e., 
cj) => 9 and 9 =>■ -itp. If there exist 6x, 62, S3, hi, h 2 such that 

ffi^i + Mi + hh 2 + 93^2 + h = -1, 

where Si, 62,63 € y, x' , y'} are sums of squares and h±,ti2 € M.[x, y, x' , y'], then 
9 = 9362 + \ < is such an interpolant for (f> and ip. In this example, applying our tool 
AiSat, we obtain in 0.025 seconds that 

In = -290.17 - 56.86^' + 1109.95a;' + 37.59y - 32.20yy' + 386.77^' + 203. 8Sy 2 + 107.91a; 2 , 
h 2 = -65.71 + Q.39y' + 244.14k' + 274. 80y + 69.33yy' - 193A2yx' - 88.18y 2 - 105.63a: 2 , 
<5i = 797.74 - 31.38^' + 466. 12y' 2 + 506.26a;' + 79. 87a;' y + 402.44a-' 2 + 104. 43j/ 

+41.09yy' - JO.liyx' + 451. 64j/ 2 + 578.94a; 2 
S 2 = 436.45, 

8 3 = 722.62 - 91.59^' + 407.17a' 2 + 69.39a' + 107.41a-'a' + 271.06a' 2 + 14.23a + 188.65aa' 

+69.33aa' 2 - 600.47aa;' - 226.01yx'y' + U2.62yx' 2 + 325.78a 2 - 156.69aV + 466.12a 2 a' 2 
+ 10.54a 2 a'a' + 595.87a 2 a;' 2 - 11.26a 3 + 41.09a 3 a' + 18.04a ; V + 451.64a 4 + 722.52a; 2 
-80.15x 2 a' + 466.12a- 2 a' 2 - 495.78a; 2 a;' + 79.87a 2 a;'a' + i02Aix 2 x' 2 + 64.57x 2 y 
+241.99a 2 a;' + 73.29x 2 yy' - 351.27a 2 aa;' + 826.70a 2 a 2 + 471.03a 4 . 

Note that <5i can be represented as 923.42(0.90 + 0.7y - O.ly' + 0.43a;') 2 + 252.84(0.42 - 
0.28i/ + 0.21y'- 0.84a;') 2 +461.69(-0.1-0.83y + 0.44y' + 0.34a;') 2 +478(-0.06 + 0.48y + 
0.87j/'+0.03a;') 2 +578.94(a;) 2 . Similarly, 6 2 and 63 can be represented as sums of squares 
also. 

Moreover, using the approach in we can prove 9 is an inductive invariant of 
the loop, therefore, error() will never be executed. 

Remark 1. Note that 9 itself cannot be generated using quantifier elimination (QE for 
short) approach in [25 1, as it contains more than thirty monomials, which means that 
there are more than thirty parameters at least in any predefined template which can be 
used to generate 9. Handling so many parameters is far beyond the capability of all the 
existing tool based on QE. However, the problem whether 9 is an inductive variant only 
contains 4 variables, therefore it can be verified using QE. The detailed comparison 
between our approach reported in this paper and QE based technique can be seen in the 
related work. 

3 Theoretical Foundations 

In this section, for self-containedness, we briefly introduce some basic notions and 
mathematical theories, based on which our approach is developed. 

4 As gi > => 52 > 0, we ignore 32 > in <j). 



Definition 1 (Interpolants). A theory T has interpolant if for all formulae cp and ip in 
the signature ofT, if (p \—-r i/l, then there exists a formula that contains only symbols 
that <p and ip share such that (p \=q- and \=j- ip. 

An interpolant of (p and —*ip is called inverse interpolant o f (p and ip, i.e., (p A 
ip |=T-L> 4> 6* an d ©Alp |=T-U where contains only the symbols that (p and ip 
share. 

Note that in practice, people like to abuse inverse interpolant as interpolant. Thus, 
as a convention, in the sequel, all interpolants are referred to inverse interpolant if not 
otherwise stated. 

Also, in what follows, we denote by x a variable vector (x\ , ■ ■ ■ , x n ) in W 1 , and by 
R[x] the polynomial ring with real coefficients in variables x. 

3.1 Problem Description 

Here, we describe the problem we consider in this paper. Let 

kt si 

Ti t = /\/y(x)t>0 and T 2l = /\ <?y(x) >' 0, (1) 

3=0 3=0 

be two semi-algebraic systems (SASs), where and gij are polynomials in R[x], 
and >ij,>ij € {=, ^, >}■ Clearly, any polynomial formula <p can be represented as a 
DNF, i.e. the disjunction of a several SASs. Let 71 = Vt=i Tit, 7a = V"=i 7~2i be two 
polynomial formulas and 71 A 7j |= JL, i.e., 71 and T2 do not share any real solutions. 
Then, the problem to be considered in this paper is how to find another polynomial 
formula / such that T\\= I and I l\Ti (= -L 

It is easy to show that if, for each t and I, there is an interpolant I t i for 7it and 722, 
then I — V"=i A"=i I tl ' s an interpolant of 71 and 72- Thus, we only need to consider 
how to construct interpolants for two SASs of the form ([1} in the rest of this paper. 

3.2 Common variables 

In the above problem description, we assume 71 and T2 share a set of variables. 
But in practice, it is possible that they have different variables. Suppose V(7i) for the 
set of variables that indeed occur in %, for i = 1,2. For each v E V(71) — V(72), if 
v is a local variable introduced in the respective program, we always have an equation 
v = h corresponding to the assignment to v (possibly the composition of a sequence of 
assignments to v); otherwise, v is a global variable, but only occurring in 71, for this 
case, we introduce an equation v = v to 75; Symmetrically, each v £ V(72) — V(71) 
can be coped with similarly. 

In the following, we show how to derive the equation v = h from the given pro- 
grams by case analysis. 

- If the given program has no recursion nor loops, we can find out the dependency 
between the variables in V(71) fl V(7a) an d the variables in V(7~j) — V(73-j) 
according to the order of assignments in the program segment, where j = 1,2. 



Clearly, we can always represent each variable in V(7j) — V(73-j) by an expres- 
sion of V(71) n V(72)- Obviously, if all expressions in the program segment are 
polynomial, the resulted expressions are polynomial either. 
- If the given program contains loops or recursion, it will become more complicated. 
So, we have to unwind the loop and represent each variable in V(7}) — V(73-j) 
by an expression of V(71) n V(72) and the number i of the iterations of the loops 
or recursions. However, the resulted expressions may not be polynomial any more. 
But as proved in ll26l . if assignment mappings of the loops in the program segment 
are solvable, the resulted expressions are still polynomial. 

Definition 2 (Solvable mapping [261). Let g € Q[x] m be a polynomial mapping, g is 
solvable if there exists a partition ofx. into subvectors of variables, X = W\ U • • • UWk, 
Wi f~)Wj — $ifi^ j, such that\/j : 1 < j < k we have 

g W] (x) = Mjwj +Pj(wi, . . 

where Mj € Q\ w i\ x \ w i\ is a matrix andPj is a vector of \wj \ polynomials in the ring 
Q[wi, . . . ,Wj—i]. For j = 1, Pi must be a constant vector, implying that g Wl isanaffine 
mapping. 



1 assume ( a+b>0) ; 1 int z = 0;int w=0; 

2 int x=a;int y=a; 2 while (w<z + 50) { 

3 while (x+y<20){ 3 a=a + l;b=b/2; 

4 x=x*a; 4 w=w+a ; z=w*b ; 

5 y=x+y*b; 5 } 

6 } 6 a s s ert (x+y+z+w>0) ; 

Code 1.2: demo example 



For example, in the Code 1.2 a, b are common variables, x, y, z, w are local vari- 
ables. Let 71 be related to the left and T2 to the right of Code 1.2 71 uses an order 
y>~x>~b>~aovi variables, and T2 uses an order z ^ w ^ d ^ a on vari- 
ables. Obviously, in every iteration of the loops, variable with higher precedence can 
only be assigned with a polynomial of variables with lower precedence. In order to 
prove the assert, we unwind the first loop i times, and obtain the values of a, 6, x, y are 
a, b, a l+1 , X)j=o +1 respectively. Similarly, unwind the second loop j times, 

and obtain the values of a, b, z, w are a + j, ^j,ja+ ii^M ; (j a _|_ U+ 1 )j respec- 
tively. Thus, in the first loop, the local variables x, y are represented by expressions of 
a, b, so are z, w in the second loop. Using such replacements, we can obtain an inter- 
polant Iij only concerning the common variables a, b w.r.t. the z-th unwinding of the 
first loop and the j-th unwinding of the second loop. Whenever we can prove that Jy is 
an invariant of Code 1 .2 then the assert is guaranteed. This is a procedure of BMC. 

In what follows, we use subvar iable to denote the above procedure to transform 
two SASs that may not share same variables to two SASs that share same variables. 



3.3 Real Algebraic Geometry 



In this subsection, we introduce some basic notions and results on real algebraic 
geometry, that will be used later. 

Definition 3 (ideal). Let X be an ideal in K[x], that is, X is an additive subgroup of 
R[x] satisfying fg G X whenever / G X and g G R[x]. Given hi, ... , h m G K[x], 

(hi, ■ ■ ■ , h m ) = 1 2j=i u jhj I Ui,.. . ,u m G K[x]| denotes the ideal generated by 

hi, ... , h m . 

Definition 4 (multiplicative monoid). Given a polynomial set P, let Mult(P) be the 
multiplicative monoid generated by P, i.e., the set of finite products of the elements of 
P (including the empty product which is defined to be 1). 

Definition5 (Cone). A cone C ofW.[x] is a subset o/E[x] satisfying the following 
conditions: (i) pi,P2 G C => p\ + P2 G C; (ii) P\,P2 G C => P1P2 G C; (Hi) p G 
M[x] p 2 G C. 

Given a set P C R[x], let C(P) be the smallest cone of M[x] that contains P. It is 
easy to see that C(0) corresponds to the polynomials that can be represented as a sum of 
squares, and is the smallest cone in M[x], i.e., { J2t=i Pi I Pi) ■ ■ • iPs € ^[ x ] }> denoted 
by SOS. For a finite set P C M[x], C(P) can be represented as: 

r 

C(P) = {q + Y^<hPi I q,qi,---,q r eC(®), Pl ,...,p r £Mult(P)}. 

i=l 

Positivstellensatz Theorem, due to Stengle [21], is an important theorem in real 
algebraic geometry. It states that, for a given SAS, either the system has a solution in 
M. n , or there exists a certain polynomial identity which bears witness to indicate that the 
system has no solutions. 

Theorem 1 (Positivestellensatz Theorem, [21|). Let (f 3 ) s J=1 , Q7fc)£=i> (M"=i defi- 
nite families of polynomials in R[x], Denote by C the cone generated by (fj)j-i, Mult 
the multiplicative monoid generated by (pfe)^. =1 , andX the ideal generated by (hi)f =1 . 
Then the following two statements are equivalent: 

( /i(x)>0, ••• , / s (x) >0, 

1. the SAS I <7i(x) ^ 0, 5t(x) ^ 0, has no real solutions; 

{ fti(xj=0, ••• , /i„(x)=0 

2. f/jere CTiif / G C, g G MwZf, h EX such that f + g 2 + h = 0. 

3.4 Semidefinite Programming 

In iBTI . Stengle did not provide a constructive proof to Theorem [T] However, Par- 
rilo in [22 23 1 provided a constructive way to obtain the witness, which is based on 
semidefinite programming. Parrilo's result will be the starting point of our method, so 
we briefly review semidefinite programming below. We use Sym n to denote the set of 
n x n real symmetric matrices, and deg(f) the highest total degree of / for a given 
polynomial / in the sequel. 



Definition 6 (Positive semidefinite matrices). A matrix M £ Sym n is called positive 
semidefinite, denoted by M >z 0, ifx T Mx > Ofor all x £ M. n . 

Definition 7 (Inner product). The inner product of two matrices A — (ay),.B = 
(bij) £ R nxn , denoted by (A, B), is defined by Tr(A T B) = E"j=i a ij b ij- 

Definition 8 (Semidefinite programming (SDP)). The standard (primal) and dual 
forms of a SDP are respectively given in the following: 

p* = inf (C, X) s.t. X t 0, (Aj,X) = bj (j = 1, . . . , m) (2) 

X GSy7n n 

rn 

d* = sup b T y s.t V VjAj + S = C, S h 0, (3) 

ye« m j=1 

where C,Ai,..., A m , S £ Sym n and b £ W n . 

There are many efficient algorithms to solve SDP such as interior-point method. We 
present a basic path-following algorithm to solve Q in the following. 

Definition 9 (Interior point for SDP). 

intF p = {X : (A u X) = bi (i = 1, . . . , m), X y 0} , 
intF d = |(y,5) : S = G ~J2 A ^ y °| ' 
intF = intFp x intFd. 

Obviously, (C,X) - b T y = (X, S) > for all (X,y, S) £ intF. Especially, we 
have d* < p* . So the soul of interior-point method to compute p* is to reduce (X, S) 
incessantly and meanwhile guarantee (X, y, S) £ intF. 



Algorithm 1: Interior_Point_Method 

input : C, Aj,bj (J = 1, . . . , m) as in ^ and a threshold c 
output: p* 

1 Given a (X, y, S) £ intF and XS = /il; 

/* /i is a positive constant and / is the identity matrix. 
*/ 

2 while fi > c do 

3 p. — 7/i; 

/* 7 is a fixed positive constant less than one */ 

4 use Newton iteration to solve (X, y, S) £ intF with XS = pi; 
s end 



Algorithm 2: Certificate-Generation 



input : {fi,...,f n },g,{hi,...,h u },b 

output: either {po, . . . ,p n } and {qi, . . . , q u } such that 

1 +Po H + Pnfn + g + qihi H + g u 7i u = 0, or NULL 

1 Let qn, (?i2, 1721, 922, • • • , q u i,qu2 G SOS with deg(gii) < b and deg(g l2 ) < febe 
undetermined SOS polynomials; 

2 Letpi, . . . ,p n £ SOS with deg(pi) < b be undetermined SOS polynomials; 

3 Let / = 1 + po + P1/1 H hPn/n + g + (qu - qi2)hi H h (<j u i - q U 2)h u ; 

4 for every monomial m £ / do 

5 Let (Q m ,Q) = coeff(/, m); 

/* Applying Lemma [lj */ 
/* coeff(/, m) the coefficient of monomial m in polynomial 

/ */ 
/* Z is a monomial vector that contains all monomials 
with coefficient 1 and degree less than or equal to 
b/2 */ 
/* po = Z T Q Z, Pl = Z T QiZ, . . . ,p n = Z T Q„Z */ 
/* Qii = Z T Qa^, 9i2 = Z T Q i2 Z, i = 1, . . . ,u */ 
/* Q = diag(l, Q , Q\, . . . ,Q n , 1, Qn,Qi2, ■ ■ ■ ,<5«i,<3«2) */ 

6 end 

7 Applying SDP software CSDP to solve whether there exists a semi-definite symmetric 
matrix Q s.t. (Q m , Q) = for every monomial m 6 / 

8 if ?Ae return o/CSDP is feasible then 

/ * qi = qu — qa * / 

9 return {p , ... ,p„} , {<?i, ... ,q u } 

10 else 

11 return NULL 

12 end 



3.5 Constructive Proof of Theorem [T] Using SDP 

Given a polynomial /(x) of degree no more than 2d, f can be rewritten as / = 
Z T QZ where Z is a vector consists of all monomials of degrees no more than d, e.g., 



[1. 



xi,x 2 , ■ ■ ■ 7 x n ,x 1 X2,X2X 3 , , and Q 



( a x 

\ 2 



2 

a„2 



2 2 ~ / 

a symmetric matrix. Note that here Q is not unique in general. Moreover, / £ C(0) 
iff there is a positive semidefinite constant matrix Q such that /(x) = Z T QZ. The 
following lemma is an obvious fact on how to use the above notations to express the 
polynomial multiplication. 

Lemma 1. Forgiven polynomials fi, ...,/„, gi, g n , assume ^™=i /».9j = E;=i c i m i 
where Q6l andrriiS are monomials. Suppose gi — Z T Q2iZ and Q2 — diag(Q2i, ■ ■ ■ , Q2 
Then there exist symmetric matrices Qu, ■ ■ ■ , Qi s such that Ci — (Qu, Q2), SILi fi9i 
J2i=i (Qui Q2) mi, in which Qu can be constructed from the coefficients of fi ,...,/„. 






Example 1. Let / = a 2 $x\ + a\\X\x 2 + a,Q 2 x 2 an d 9 = ^00 + frio^i + boiX 2 . Then, 
f9=(Qu, Q2) xj + {Qi2, Q2) x 1 x 2 + {Qi3, Q 2 ) xj+(QiA, Q2) x\x\ (Q15, Q2) x\x 2 + 
(Qie, Q2) x\ + {Qn, Q 2 ) x\ , where 

(a 02 0\ 

O2 = I T h On = I I Qi2 = I I , Q13 = 

\ 00/ 

(Q £02 ( 
5f C 
c 

Back to Theorem [TJ We show how to find / G C, g G Mult, h £ 1 such that 
/ + g 2 + h = via SDP solving. First, since / G C, f can be written as a sum 
of the products of some known polynomials and some unknown SOSs. Second, h E 
l({hi, . . . , h u }) is equivalent to h = hipi + ■ ■ ■ + h u p u , which is further equivalent 

to h = hi(qu - q 12 ) H h h u (q ul - q u2 ), where p^q^ G R[x] and q l} E SO^] 

Third, fix an integer d > 0, let g = (nf =1 gi) d , and then / + g 2 + h = can be written 
as Yli=i fi^u where I is a constant integer, f- & R[x] are known polynomials and Si G 
SOS are undermined SOS polynomials. Therefore, Theorem[TJis reduced to fixing a 
sufficiently large integer d and finding undetermined SOS polynomials S{ occurring in 
/, h with degrees less than or equal to deg(g 2 ), which satisfies f+g 2 + h = 0. Based on 
Lemma[T[ this is a SDP problem of form Q. The constraints of the SDP are of the form 
(Aj,X) = 0, where A, and X correspond to Qij and Q 2 in Lemma [T| respectively. 
And Q 2 is a block diag matrix whose blocks correspond to the undetermined SOS 
polynomials in the above discussion. That is, 

Theorem 2 (|22|). Consider a system of polynomial equalities and inequalities of the 
form in Theorem [7] Then the search for bounded degree Positivstellensatz refutations 
can be done using semidefinite programming. If the degree bound is chosen to be large 
enough, then the SDPi will be feasible, and the certificates can be obtained from its 
solution. 

Algorithm [2] is an implementation of Theorem [2] and we will invoke Algorithm [2] 
as a subroutine later. Note that Algorithm|2]is a little different from the original one in 
, as here we require that / has 1 as a summand for our specific purpose. 



4 Synthesizing Non-linear Interpolants in General Case 

As discussed before, we only need to consider how to synthesize interpolants for 
the following two specific SASs 

r/ 1 (x)>o,...,/ sl ( x )>o ! r / sl+1 (x)>o,...,/ s (x)>o, 

Tx = I .91 (x) ^0,...,. 9tl (x) ^0, T 2 = I 5 t 1+ i(x)^0,..., fft (x)^0, (4) 
[ hi(x) = 0, . .., h Ul (x) =0 [ h Ul+ i(x) = 0, . . .,h u (x) = 

where 71 and T 2 do not share any real solutions. 

5 For example, let q a = (\pi + l) 2 ,<?i2 = ( jPi - l) 2 - 



By Theorems [jM] there exist / e C({/i, . . . , / s }), g e Mult({gi, . . .,g t }) and 
/i e ,h u }) such that / + g 2 + h = 0, where 



h = q 1 h 1 -\ h q ul h Ul H h g„/i u , 

/ = Po+Pifi H hp s / s +P12/1/2 H hpi... s /i . . . / s . 

in which and are in SOS. 

If / can be represented by three parts: the first part is an SOS polynomial that 
is greater than 0, the second part is from C({f%, . . . , f Sl }), and the last part is from 

■ • ■ : fs}), i-e., / = Po+J2vC{i,..., Sl } P u(^iS«/i)"l"X/uC{si + l s} Pv{ni£ V fi), 

where Vx e R n .p (x) > and p v £ SOS. Then let 

/Ti = X! Pv n i£vfi, h Tl = qihi H V q Ul h Ul , 

dC1,,,..si 

/r 2 = 51 Pvllitvfi, h T2 = h~h Tll 

q = fn+ g 2 + h Tl + -~ = -(/r 2 + V 2 ) - 

Obviously, we have VxG 71.(7(x) > and Vx£72-g(x) < 0. Thus, let 7 = </(x) > 0. 
We have Ti |= I and 7 A 72 |=±. 

Notice that because the requirement on / cannot be guaranteed in general, the above 
approach is not complete generally. We will discuss under which condition the require- 
ment can be guaranteed in the next section. We implement the above method for syn- 
thesizing non-linear interpolants in general case by Algorithm[3] 

Example 2. Consider 

( x\ + x\ + x\ - 2 > 0, ( -Zx\ - Ax\ - 10x2 + 20 > 0, 

T\ = < xi + x 2 + x 3 7^ 0, and T2 = I 2x 1 + 3x 2 - 4x 3 7^ 0, 

[ 1.2x1 + x 2 + x i x 3 = { xl + x% — x 3 — 1 = 




Algorithm 3: SN_Interpolants 



input : 71 and T2 of the form Q, b 
output: An interpolant I or NULL 

1 g := nl =x gl 

I h 

2 g \= g de s(s) J 

3 {/ii} : = { n ievfi foru C {1, .. . ,si}} ; 

4 {/i 2 } : = {n iev fi foru C {si + l,...,s}}; 
s Vi = V({f tl }U{hi,...,h Ul }); 

/* Get all variables in polynomials system */ 

6 v 2 = v({f t2 }u{h ui+1 ,...,h u })- 

7 V = V! n V 2 ; 

8 ({/ti>/t 2 }iffi {^i. • ■ • > fe u}) = subvariable( {/ tl , f t2 }, g, 
{hi, . .., h u },V,{hi, .. ., h u })\ 

/* Replace every uncommon variable v by a polynomial h where 
v — h as described in Subsection 3.2 */ 

9 sdp:=Certif icate_Generat ion({/t 1 , ft 2 }, g, {hi, ■ ■ ■ , h u } , 6) 

10 if sdp = NULL then 

11 return NULL 

12 else 

13 I := I + Z„c{i,..., sl }P^e«/; + <7i^i H h g Ul /iui + g > 0; 

14 return 7; 

15 end 



Clearly, 71 and T2 do not share any real solutions, see Fig. [I] (left) By set- 
ting 6 = 2, after calling Certificate-Generation, we obtain an interpolant I 
with 30 monomials -14629.26 + 2983. 44x 3 + 10972. 97x1 + 297.62x2 + 297.64x 2 x 3 + 
0.02x22;! + 9625.61a;| - 1161.80x|x 3 + 0.01x|x! + 811.93xi + 2745. Mx^ - 10648. llxi + 
3101.42xix 3 + 8646. 17xix| + 511. 84xiX2-1034. 31xix 2 x 3 + 0.02xix 2 xi + 9233. 66xix| + 
1342.55x1x1x3 - 138.70x!x| + 11476.61xf - 3737.70x?x 3 +4071.65x?x§ - 2153.00x?x 2 + 
373.14x?x 2 x 3 + 7616. 18x?x| + 8950.77x? + 1937.92xf x 3 - 64.07x?x 2 +4827.25xf , whose 
figure is depicted in Fig.[T](right). □ 



5 A Complete Algorithm Under Archimedean Condition 

Our approach to synthesizing non-linear interpolants presented in Section [4] is in- 
complete generally as it requires that the polynomial / in C({/i, . . . , / s }) produced by 
Algorithm [2] can be represented by the sum of three polynomials, one of which is posi- 
tive, the other two polynomials are respectively from C({/i, . . . ,/ Sl }) an dC({/ Sl +i, ■ ■ ■ , 
In this section, we show, under Archimedean condition, the requirement can be indeed 
guaranteed. Thus, our approach will become complete. In particular, we shall argue 
Archimedean condition is a necessary and reasonable restriction in practice. 



6 For simplicity, we do not draw xi + xi + X3 7^ 0, nor 2xi + 3x2 — 4x3 7^ in the figure. 



5.1 Archimedean Condition 

To the end, we need more knowledge of real algebraic geometry. 
Definition 10 (quadratic module). For g\, ... , g m G E[x], the set 

m 

M(gi, ■■■,g m ) = {5 + Yl 6 M I S °> § j e C W> < 5 > 

3=1 

is called the quadratic module generated by g\ , . . . , g m . A quadratic module Ai is 
called proper if —1 ^ Ai (i.e. Ai ^ K[x]J. A quadratic module Ai is maximal 
for any p £ l[x] fl Ai, Ai U {p} is not a quadratic module. 

In what follows, we will use — Ai to denote {— p | p G Ai} for any given quadratic 
module Ai . 

The following results are adapted from [27 | and will be used later, whose proofs 
can be found in ll27l . 

Lemma 2 ( Il2ll27l ). 

1) If Ai C K[x] is a quadratic module, then I = Ai fl — Ai is an ideal. 

2) If M C R[x] a maximal proper quadratic module, then Ai U — .M = R[x]. 

3) {x G R™ | /(x) > 0} is a compact se^for some f G M({fi, . . . , /«}) ;7J 

Vp G R[x], 3n G N.n ±p G . . . , /.). (6) 

Definition 11 (Archimedean). For gi,...,g m G R[x], f/ie quadratic module Ai (g±, . . . , 
is 50/0" fo foe Archimedean i/f/ze condition (J6| holds. 

Let 

Ti = A(x) > 0, . . . , / Sl (x) > and T 2 = / S1+1 (x) > 0, . . . , / s (x) > (7) 

be two SASs, where {/j(x) | i = 1, . . . , s} contains constraints on the upper and lower 
bounds of every variable X4, and 71 and T2 do not share real solutions. 

Remark 2. Regarding {/1, . . . , f s } in Q, as every variable is bounded, assume N — 
Y^i—i x i £ {fl: ■ ■ ■ 1 fs} for a const N, then Ai(f\, f s ) is Archimedean. 

Lemma 3. [21,27)1 Let Ai C R[x] foe a maximal proper quadratic module which is 
Archimedean, J = Ai fl — A4, anc/ / G R[x], f/ien f/iere exists a 6 R iMcfo f/iflf 
/-AG/. 

Lemma 4. 7f7 is an ideal and there exists el = (ax , . . . , a„) G R n iMcfo f/za? X, — a, G 7 
/or i = 1, . . . , n, then for any / G R[x] , / — /(a) G 7. 

Proof. Because x$ — 0,4 G J for i = 1, . . . , n, (xi — a±, . . . ,x n — a n ) C 7. For any 
/ G R[x], (xi — ai, . . . ,x n — a n ) is a radical ideaQand (/ — /(a))(a) = 0, so / — 
/(a) G (x x - a\, . . . ,x n - On) C 7. □ 

7 5 is a compact set in R n iff S 1 is a bounded closed set. 

8 Ideal 7 is a radical ideal if 7 = yl — {/|/ G 7 for some integer > 0}. 



Theorem 3. Suppose {/i(x), . . . , / s (x)} is given in |7|). If > 0) is unsatisfi- 

able, then -1 e M(fi, ■ ■ ■ , / s ). 

Proof. By Remark[2] .M (/i , . . . , / s ) is Archimedean. Thus, we only need to prove that 
the quadratic module -M(/i, . . . , f s ) is not proper. 

Assume A4(fi, ■ ■ ■ , f s ) is proper. By Zorn's lemma, we can extend A4(fi, . . . , f s ) 
to a maximal proper quadratic module Ai 2 Ad(f\, . . . , f s ). As M.(f\, ■ ■ ■ , f s ) is 
Archimedean, .M is also Archimedean. By Lemma[3] there exists a = (ai, . . . , a n ) £ 
K" such that Xi-a 4 e / = MO-M for all i G {1, . . . , n}. FromLemma|4] /-/(a) £ 
7 for any / e M[x]. In particular, for f — fj, we have /*/ (a) = /j — (/j — fj(a)) £ Ai 
since / 3 £ M(fi, ■ ■ ■ , f s ) Q M and — (fj - fj(a)) £ M, which implies fj(a) > 0, 
for j = 1, . . . , s. This contradicts to the unsatisfiability of A<=i (fi — 0)- ^ 

By Theorem|3]we have —1 £ M.(fi, . . . , / s ). So, there exist oo, ■ ■ ■ , cr s £ C(0) 
such that -1 = cr + a\f\ -\ h <r Sl f Sl + a Sl+1 f Sl+1 H h / S <7 S . It follows 

- + ^si + l/i'i + l H 1" <?sfs) = 2 + CT o + cx/i H h cr Sl / Sl . (8) 

Let q(x) = I + cr + 01/1 + • • • + c Sl / Sl , we have Vx e 71.g(x) > and Vx e 
72-9(x) < 0. Thus, let 7 = g(x) > 0. According to Definition[T] 7 is an interpolant of 
71 and 72- So, under Archimedean condition, we can revise Algorithm[3]as Algorithm|4] 



Algorithm 4: RSN_Interpolants 
input : Ti and T2 as in (J7J, 



,h u } 



/* h ui +i,...,h u are the equality occur in T2 
output: 7 

1 b=2; 

2 Vi = V({/i,...,/^}); 

/* Get all variables of 7i 

3 v 2 = v({f si+u ...,f u }); 

4 V = Vi(lV 2 ; 

5 {/si+i>- • • ,/»})=subvariable({/ Sl +i, . . .,/«}),V, {Vli- ■ • , K}); 

/* Replaceing every uncommon variable v by polynomial h 
where v>h and v < ft as described in Section 3.2 

6 while true do 
sdp=Certif icate_Generation({/i, . . . , / s },0,{},6); 
if sdp / NULL then 

I={h + ULiPifi>0h 
/'=subvariable(7,V, . . . , 

return 7'; 

else 

I b=b+2; 
end 



7 




8 
9 




10 




11 




12 




13 




14 




15 end 



*/ 



*/ 



*/ 



Example 3. Let & = /\? =1 Xi > —2 A > -2, /i = — arf - 4^2 - + 2, 
/a = x\ — x\—x\x^ — \, f 3 = — x\ — ^x\ — x|+3xi^2+0.2, / 4 = — a;f +x|+x 1 a;3+l. 
Consider 71 = # A f x > A f 2 > and 7~ 2 = & A / 3 > A f 4 > 0. Obviously, 71 A T 2 
is unsatisfiable, see Fig. [2] (left). 




Fig. 2: Example[3] 

By applying RSN_Interpolants, we can get an interpolant as — 33.7255a-i + 
61.1309xfx 2 + 4.6818x?a;3 - 57.927a;?a;l + 13.4887a;fa;2a;3 - 48. 9983a;? zi - 8.144a;? - 
48.1049o;ia;|-6.7143a;ia;la;3+29.8951a;iX2a;l+61.5932a;ia;2+0.051659a;ix|-0.88593a;iX3- 
34.7211xf-7.8128a;|a;3-71.9085a;|a;i-60. 5361i-|-l. 6845a;2a;|-0.5856a; 2 a;3-15.2929a;f- 
9.75633;! + 6.7326, which is depicted in Fig [5] (right). In this example, the final value of 
6 is 2. □ 



5.2 Discussions 

1. Reasonability of Archimedean condition: Considering only bounded numbers can 
be represented in computer, so it is reasonable to constraint each variable with upper and 
lower bounds in practice. Not allowing strict inequalities indeed reduce the expressive- 
ness from a theoretical point of view. However, as only numbers with finite precision 
can be represented in computer, we always can relax a strict inequality to an equiva- 
lent non-strict inequality in practice. In a word, we believe Archimedean condition is 
reasonable in practice. 

2. Necessity of Archimedean condition: In Theorem [5] Archimedean condition is 
necessary. For example, let 71 = {x% > 0, x 2 > 0} and T2 = {— X1X2 — 1 > 0}. 
Obviously, 71 A T2 = is not Archimedean and unsatisfiable, but 

Theorem 4. — 1 ^ M(xi,x 2 , — Xxx 2 — 1). 

Proof. Suppose these exist 60,61,62,63 £ C(0) such that h = 60 + x\6\ + X262 — 
{xix 2 +l)6 3 = -1. Let cose?" °xf °, Cxxf 1+1 xf\ c 2 x\ a * xf 2+1 , and c z x\ a3+1 xf 3+1 
be the leading terms of 8$, x\6\, X262 and {x\x 2 + 1)#3, respectively, according to the 
total degree order of monomials, where c, > and at , 6; £ N. Obviously, the four 
terms are pairwise different. So, the leading term of h must be one of the four terms if 
they are not zero. This, together with h = —1, imply that c\ — c 2 = C3 = and thus 
8\ = 62 = 63 = 0. Therefore, 6q = —1, a contradiction. □ 



6 Correctness and Complexity Analysis 



The correctness of the algorithm SN_Interpolants is obvious according to 
Theorem [2] and the discussion of Section 4. Its complexity just corresponds to one 
iteration of the algorithm RSN_Interpolants. The correctness of the algorithm 
RSN_Interpolants is guaranteed by Theorem[2]and Theorem[3] The cost of each 
iteration of RSN_Interpolants depends on the number of the variables n, the 
number of polynomial constraints u, and the current value of bf. The size of X in 
^ is tt(" + ^ ^ 2 ) and the m in is (™^ t b/ )- So, the complexity of applying interior 
method to solve the SDP is polynomial in u( n+ ^^ 2 ) ( n+b/ ). Hence, the cost of each 
iteration of RSN.Interpolants is u( n+b n f/2 ) ( n+ J' f ). Therefore, the total cost of 
RSN.Interpolants is &yu( n+ ^/ 2 ) ( n+ J' f )- For a given problem, n,u are fixed, so 
the complexity of the algorithm becomes polynomial in bf. The complexity of Algo- 
rithm SN_Interpolants is the same as above discussions, except that the number 
of polynomial constraints is about 2 Sl + 2 s ~ Sl . 

As indicated in (23|, there are upper bounds on bf, which are at least triply expo- 
nential. So our approach can enumerate all possible instances, but can not be done in 
polynomial time. 



7 Implementation and Experimental Results 

We have implemented a prototypical tool of the algorithms described in this pa- 
per, called AiSat, which contains 6000 lines of C++ codes. AiSat calls Singular 
1 28 1 to deal with polynomial input and CSDP to solve SDPs. In AiSat, we design a 
specific algorithm to transform polynomial constraints to matrices constraints, which 
indeed improves the efficiency of our tool very much, indicated by the comparison with 
SOSTOOLS l|29l (see the table below). As a future work, we plan to implement a new 
SDP solver with more stability and convergence efficiency on solving SDPs. 

In the following, we report some experimental results by applying AiSat to some 
benchmarks. 



The first example is from 1301 . see the source code in Code 1.3 We show its cor- 



rectness by applying AiSat to the following two possible executions. 

- Subproblem 1: Suppose there is an execution starting from a state satisfying the 
assertion at line 13 (obviously, the initial state satisfies the assertion), after — > 6 — > 
7—5-8^9^11^12^13, ending at a state that does not satisfy the assertion. 
Then the interpolant synthesized by our approach is 716. 77+1326. 74(ya)+1.33(ya) 2 - 
433.90(?/a) 3 + 668.16(xa)-155.86(a;a)(2/a) + 317.29(xa)(ya) 2 + 222.00(a;a) 2 + 
592.39(a;a) 2 (?/a) + 271.11(xa) 3 , which guarantees that this execution is infeasible. 

- Subproblem 2 : Assume there is an execution starting from a state satisfying the 
assertion at line 13, after — > 6 — » 7 — » 8 — > 10 — > 11 — > 12 — > 13, ending at a state 
that does not satisfy the assertion. 

The interpolant generated by our approach is 716.95 + 1330.91(ya) + 67.78(ya) 2 + 
551.51(ya) 3 + 660.66(:ra) - 255.52{xa)(ya) + 199.84(xa)(ya) 2 + 155.63(:m) 2 + 386.87 
(xa) 2 (ya) + 212.41(a;a) 3 , which guarantees this execution is infeasible either. 



1 


int main () { 


1 


vc :=0; 


2 


int x , y ; 


2 


/* the initial veclocity */ 


3 


int xa := 0; 


3 


fr :=1000; 


4 


int ya := 0; 


4 


/* the initial force */ 


c 

3 


while (nondet()) { 


c 

3 


ac : =U . 0005 * tr ; 


6 


x := xa + 2*ya; 


6 


/* the initial acceleration */ 


7 


y := — 2*xa + ya ; 


7 


while ( 1 ) { 


8 


x + + ; 


8 


fa:=0.5418*vc*vc; 


9 


if (nondet()) y= y+x ; 


9 


/* the force control */ 


10 


else y : = y— x ; 


10 


fr :=1000-fa ; 


11 


xa := x — 2*y ; 


11 


ac : = 0.0005* fr ; 


12 


ya := 2*x + y;} 


12 


vc : = vc + ac ; 


13 


assert (xa + 2*ya >= 0) ; 


13 


assert (vc <49.61) ; 


14 


return 0; 


14 


/* the safety velocity */ 


15 } 




15 


} 




Code 1.3: ex 1 




Code 1 .4: An accelerating car 



The second example accelerate (see Code 1.4 1 is from [20|. Taking the air 
resistance into account, the relation between the car's velocity and the physical drag 
contains quadratic functions. Due to air resistance the velocity of the car cannot be 
beyond 49.61to/s, which is a safety property. Assume that there is an execution (vc < 
49.61) -> 8 -> 10 -> 11 -> 12 -> 13(uc > 49.61). By applying Applying AiSat, we 
can obtain an interpolant — 1.3983uc + 69.358 > 0, which guarantees vc < 49.61. So, 
accelerate is correct, we can synthesize an interpolant, which guarantees the safety 
property. 

The last example logistic is also from |20|. Mathematically, the logistic loop 
is written as x n+ i = rx n (l — x n ), where < x n < 1. When r = 3.2, the logistic 
loop oscillates between two values. The verification obligation is to guarantee that it is 
within the safe region (0.79 < iAi < 0.81) V (0.49 < x A x < 0.51). By applying 
AiSat to the following four possible executions, the correctness is obtained. 



- Subproblem 1: {x > 0.79 A x < 0.81} logistic {x > 0.51} is invalidated by 
the synthesized interpolant 108.92 - 214.56x > 0. 

- Subproblem 2: {x > 0.79 Ax < 0.81} logistic {x < 0.49} is outlawed by the 
synthesized interpolant -349.86 + 712. 97x > 0. 

- Subproblem 3: {x > 0.49 A x < 0.51} logistic {x > 0.81} is excluded by the 
generated interpolant 177.21 - 219.40a; > 0. 

- Subproblem 4: {x > 0.49 A x < 0.51} logistic {x < 0.79} is denied by the 
generated interpolant -244.85 + 309.31a; > 0. 



Some experimental results of applying AiSat to the above three examples on a 
desktop (64-bit Intel(R) Core(TM) i5 CPU 650 @ 3.20GHz, 4GB RAM memory and 
Ubuntu 12.04 GNU/Linux) are listed in the table below. Meanwhile, as a comparison, 
we apply the SOSTOOLS to the three examples with the same computer. 



Benchmark 


#Subporblems 


AiSat (milliseconds) 


S0ST00LS (milliseconds) 


exl 


2 


60 


3229 


accelerate 


1 


940 


879 


logistic 


4 


20 


761 



8 Related work 

In Introduction, we have introduced many work related to interpolant generation and 
its application to program verification. In this section, we will mention some existing 
work on program verification, so that we give a comparison between our approach and 
them. 

Work on program verification can date back to the late sixties (or early seven- 
ties) of the 20th century when the so-called Floyd-Hoare-Naur's inductive assertion 
method i3 11321331 was invented, which was thought as the dominant approach on auto- 
matic program verification. The method is based on Hoare Logic [32], by using pre- and 
post- conditions, loop invariants and termination analysis through ranking functions, 
etc. Therefore, the discovery of loop invariants and ranking functions plays a central 
role in proving the correctness of programs and is also thought of as the most challeng- 
ing part of the approach. 

Since then, there have been lots of attempts to handle invariant generation of pro- 
grams, e.g. [34 35 36 37 1, but only with a limited success. Recently, due to the advance 
of computer algebra, several methods based on symbolic computation have been ap- 
plied successfully to invariant generation, for example the techniques based on abstract 
interpretation [38 39 40 41 ], quantifier elimination [42 43 25 1 and polynomial algebra 
H44I45I26I46I . 

The basic idea behind the approaches based on abstract interpretation is to per- 
form an approximate symbolic execution of a program until an assertion is reached that 
remain unchanged by further executions of the program. However, in order to guaran- 
tee termination, the method introduces imprecision by use of an extrapolation operator 
called widening/narrowing. This operator often causes the technique to produce weak 
invariants. Moreover, proposing widening/narrowing operators with certain concerns of 
completeness is not easy and becomes a key challenge for abstract interpretation based 
techniques B8I391 . 

In contrast, approaches by exploiting the theory of polynomial algebra to discover 
invariants of polynomial programs were proposed in [44 45 26 46 1. In [44 1, Mueller- 
Olm and Seidl applied the technique of linear algebra to generate polynomial equa- 
tions of bounded degree as invariants of programs with affine assignments. In B45I26L 
Rodrigez-Carbonell and Kapur first proved that the set of polynomials serving as loop 
invariants has the algebraic structure of ideal, then proposed an invariant generation al- 
gorithm by using fixpoint computation, and finally implemented the algorithm by the 
Grobner bases and the elimination theory. The approach is theoretically sound and com- 
plete in the sense that if there is an invariant of the loop that can be expressed as a con- 
junction of polynomial equations, applying the approach can indeed generate it. While 
in ]46ll . the authors presented a similar approach to finding invariants represented by 
a polynomial equation whose form is priori determined (called templates) by using an 



extended Grobner basis algorithm. The complexity of the above approaches are double 
exponential as Grobner base technique is adopted. 

Compared with polynomial algebra based approaches that can only generate invari- 
ants represented as polynomial equations, Colon et al in |42| proposed an approach to 
generate linear inequalities as invariants for linear programs, based on Farkas' Lemma 
and nonlinear constraint solving. The complexity depends on the complexity of linear 
programming, which is in general is polynomial in the number of variables. 

In addition, Kapur in ||43ll proposed a very general approach for automatic genera- 
tion of more expressive invariants by exploiting the technique of quantifier elimination, 
and applied the approach to Presburger Arithmetic and quantifier-free theory of con- 
junctively closed polynomial equations. Theoretically speaking, the approach can also 
be applied to the theory of real closed fields, but Kapur also pointed out in B31l that this 
is impractical in reality because of the high complexity of quantifier elimination, which 
is doubly exponential flTl in the number of quantifiers. While in |25|, we improved 
Kapur's approach by using the theory of real root classification of SASs [48], with the 
complexity singly exponential in the number of variables and doubly exponential in the 
number of parameters. 

Comparing with the approaches based on polynomial algebra, or Farkas' Lemma or 
Grobnes, our approach is more powerful, also more efficient except for Farkas' Lemma 
based approach. Comparing with quantifier elimination based approach [43 25], our 
approach is much more efficient, even according to the complexity analysis of quantifier 
elimination given in ll49ll . which is doubly exponential in the number of the quantifier 
alternation, and becomes singly exponential in the number of variables and constraints 
in our setting 

9 Conclusion 

The main contributions of the paper include: 

- We give a sound but not inomplete algorithm SN_Interpolants for the gener- 
ation of interpolants for non-linear arithmetic in general. 

- If the two systems satisfy Archimedean condition, we provide a more practical 
algorithm RSN_Interpolants, which is not only sound but also complete, for 
generating Craig interpolants. 

- We implement the above algorithms as a protypical tool AiSat, and demonstrate 
our approach by applying the tool to some benchmarks. 

In the future, we will focus on how to combine non-linear arithmetic with other 
well-established decidable first order theories. In particular, we believe that we can use 
the method of [51 20| to extend our algorithm to uninterpreted functions. To investigate 
errors caused by numerical computation in SDP is quite interesting. In addition, to 
investigate the possibility to apply our results to the verification of hybrid systems is 
very significant. 

9 Hoon Hong in [ 50 1 pointed out the existing algorithms for the existential real theory which 
are singly exponential in the number of variables is far from realization, even worse than the 
general algorithm for quantifier elimination. 



References 



1. Clarke, E.M., Emerson, E.A.: Design and synthesis of synchronization skeletons using 
branching-time temporal logic. In: Logic of Programs. (1981) 52-71 

2. Nipkow, T., Wenzel, M., Paulson, L.C.: Isabelle/HOL: a proof assistant for higher-order 
logic. Springer- Verlag, Berlin, Heidelberg (2002) 

3. Owre, S., Rushby, J., Shankar, N.: Pvs: A prototype verification system. In: CADE. Volume 
607ofLNCS. (1992)748-752 

4. Cousot, P., Cousot, R.: Abstract interpretation: a unified lattice model for static analysis of 
programs by construction or approximation of fixpoints. In: Proceedings of POPL'77. (1977) 
238-252 

5. Biere, A., Cimatti, A., Clarke, E., Zhu, Y.: Symbolic model checking without bdds. In: 
TACAS'99. Volume 1579 of LNCS. (1999) 193-207 

6. Clarke, E., Grumberg, O., Jha, S., Lu, Y., Veith, H.: Counterexample-gided abstraction re- 
finement. In: CAV'00. Volume 1855 of LNCS. (2000) 154-169 

7. Nieuwenhuis, R., Oliveras, A., Tinelli, C: Solving sat and sat modulo theories: From an 
abstract davis-putnam-logemann-loveland procedure to dpll(t). J. ACM 53(6) (2006) 937- 
977 

8. Nelson, G., Oppen, D.C.: Simplification by cooperating decision procedures. ACM Trans. 
Program. Lang. Syst. 1(2) (October 1979) 245-257 

9. Davis, M., Logemann, G., Loveland, D.: A machine program for theorem-proving. Commun. 
ACM 5(7) (1962) 394-397 

10. Craig, W.: Linear reasoning: A new form of the herbrand-gentzen theorem. J. Symb. Log. 
22(3) (1957) 250-268 

1 1. de Moura, L., Bjrner, N.: Z3: An efficient smt solver. In: TACAS'08. Volume 4963 of LNCS. 
(2008) 337-340 

12. McMillan, K.L.: Interpolation and sat-based model checking. In: CAV'03. Volume 3920 of 
LNCS. (2003) 1-13 

13. McMillan, K.L.: An interpolating theorem prover. Theor. Comput. Sci. 345(1) (2005) 101- 
121 

14. Graf, S., Saidi, H.: Construction of abstract state graphs with pvs. In: CAV'97. Volume 1254 
of LNCS. (1997)72-83 

15. Henzinger, T.A., Jhala, R., Majumdar, R., McMillan, K.L.: Abstractions from proofs. In: 
POPL04. (2004) 232-244 

16. McMillan, K.L.: Lazy abstraction with interpolants. In: CAV'06. Volume 4144 of LNCS. 
(2006) 123-136 

17. Jung, Y, Lee, W., Wang, B.Y, Yi, K.: Predicate generation for learning-based quantifier-free 
loop invariant inference. In: TACAS'll. Volume 6605 of LNCS. (2011) 205-219 

18. Yorsh, G, Musuvathi, M.: A combination method for generating interpolants. In: CADE'05. 
Volume 3632 of LNCS. (2005) 353-368 

19. Rybalchenko, A., Sofronie-Stokkermans, V.: Constraint solving for interpolation. J. Symb. 
Comput. 45(11) (2010) 1212-1233 

20. Kupferschmid, S., Becker, B.: Craig interpolation in the presence of non-linear constraints. 
In: FORMATS' 11. Volume 6919 of LNCS. (2011) 

21. Bochnak, J., Coste, M., Roy, M.F: Real Algebraic Geometry. Springer (1998) 

22. Parrilo, P.A.: Structured semidefinite programs and semialgebraic geometry methods in ro- 
bustness and optimization. PhD thesis, California Inst, of Tech. (2000) 

23. Parrilo, PA.: Semidefinite programming relaxations for semialgebraic problems. Mathemat- 
ical Programming 96 (2003) 293-320 

24. Vandenberghe, L., Boyd, S.: Semidefinite programming. SIAM Review 38(1) (1996) 49-95 



25. Chen, Y., Xia, B., Yang, L., Zhan, N.: Generating polynomial invariants with discoverer and 
qepcad. In: Formal Methods and Hybrid Real-Time Systems. Volume 4700 of LNCS. (2007) 
67-82 

26. Rodriguez-Carbonell, E., Kapur, D.: Generating all polynomial invariants in simple loops. 
Journal of Symbolic Computation 42 (2007) 443^176 

27. Laurent, M.: Sums of squares, moment matrices and optimization over polynomials. In: 
Emerging Applications of Algebraic Geometry. Volume 149 of The IMA Volumes in Math- 
ematics and its Applications. (2009) 157-270 

28. Greuel, G.M., Pfister, G, Schonemann, H.: Singular: a computer algebra system for polyno- 
mial computations. ACM Commun. Comput. Algebra 42(3) (2009) 180-181 

29. Prajna, S., Papachristodoulou, A., Seiler, P., Parrilo, PA.: SOSTOOLS: Sum of squares 
optimization toolbox for MATLAB. (2004) 

30. Gulavani, B., Chakraborty, S., Nori, A., Rajamani, S.: Automatically refining abstract inter- 
pretations. In: TACAS'08. Volume 4963 of LNCS. (2008) 443-458 

3 1 . Floyd, R.W. : Assigning meanings to programs. In: Proc. Symphosia in Applied Mathematics 
19. (1967) 19-37 

32. Hoare, C: An axiomatic basis for computer programming. Comm. ACM 12(10) (1969) 
576-580 

33. Naur, P.: Proofs of algorithms by general snapshops. BIT 6 (1966) 310-316 

34. Wegbreit, B.: The synthesis of loop predicates. Communications of the ACM 17(2) (1974) 
102-112 

35. German, S., Wegbreit, B.: A synthesizer of inductive assertions. IEEE Transactions on 
Software Engineering 1(1) (1975) 68-75 

36. Katz, S., Manna, Z.: Logical analysis of programms. Communications of the ACM 19(4) 
(1976) 188-206 

37. Karr, M.: Affine relationships among variables of a program. Acta Informatica 6 (1976) 
133-151 

38. Cousot, P., Halbwachs, N.: Automatic discovery of linear restraints among the variables of 
a program. In: ACM POPL78. (1978) 84-97 

39. F. Besson, T.J., Talpin, J.P: Polyhedral analysis of synchronous languages. In: SAS'99. 
Volume 1694 of LNCS. (1999) 51-69 

40. Rodriguez-Carbonell, E., Kapur, D.: An abstract interpretation approach for automatic gen- 
eration of polynomial invariants. In: SAS'04. Volume 3148 of LNCS. (2004) 280-295 

41. Cousot, P.: Proving program invariance and termination by parametric abstraction, lan- 
grangian relaxation and semidefinite programming. In: VMCAF05. Volume 3385 of LNCS. 
(2005) 1-24 

42. M. Colon, S.S., Sipma, H.: Linear invariant generation using non-linear constraint solving. 
In: CAV'03. Volume 2725 of LNCS. (2003) 420-432 

43. Kapur, D.: Automatically generating loop invariants using quantifier elimination. In: Intl. 
Conf. on Applications of Computer Algebra (AC A 04). (2004) 

44. Miiller-Olm, M., Seidl, H.: Precise interprocedural analysis through linear algebra. In: ACM 
POPL04. (2004) 330-341 

45. Rodriguez-Carbonell, E., Kapur, D.: Automatic generation of polynomial loop invariants: 
algebraic foundations. In: ISSAC'04. (2004) 

46. S. Sankaranarayanan, H.S., Manna, Z.: Non-linear loop invariant generation using grobner 
bases. In: ACM POPL04. (2004) 318-329 

47. J. H., D., Heintz, J.: Real elimination is doubly exponential. J. of Symbolic Computation 5 
(1988) 29-37 

48. Xia, B., Yang, L.: An algorithm for isolating the real solutions of semi-algebraic systems. J. 
Symbolic Computation 34 (2002) 461-477 



49. Brown, C.W.: The Complexity of Quantifier Elimination and Cylindrical Algebraic Decom- 
position. United States Naval Academy 

50. Hong, H.: Comparison of several decision algorithms for the exis- tential theory of the reals 
(1991) 

51. Sofronie-Stokkermans, V.: Interpolation in local theory extensions. In: Automated Reason- 
ing. (2006) 235-250 



